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INTRODUCTION 


BACKGROUND 

Maintaining  structural  Integrity  of  thin  section  aircraft  components  when 
they  may  contain  potential  flaws  Is  the  primary  objective  of  damage  tolerance 
evaluations.  Because  Integrity  of  aircraft  structures  may  become  challenged 
by  severe  and  rapid  heating,  especially  when  the  resul  tl  ng  thermal  loads  are 
superimposed  with  the  aerodynamic  loads  of  normal  flight,  accurate  damage 
tolerance  analysis  must  oonsider  the  effect  of  thermal  environment  on 
structural  reliability  (1,  2). 

High  Intensity  heating  of  aircraft  structures  may  be  produced  by  aerodynamic 
heating,  by  laser  Irradiation,  or  by  localized  Intense  fire.  Structural 
reliability  will  be  reduced  by  the  development  of  thermal  stresses  as  well  as 
the  degradation  of  strength  properties  at  elevated  temperatures  caused  by  the 
heating.  The  failure  modes  resulting  from  these  conditions  Include  localized 
melting  or  complete  burn  through.  Such  conditions  could  result  In  brittle  or 
ductile  fracture,  or  buckling  of  thin  sections  depending  on  the  nature  of  the 
combined  stress  state  and  the  type,  size,  and  location  of  existing  flaws. 

Fully  plastic  failure  of  unflawed  aluminum  panels  by  plastic  collapse  has  been 
studied  to  the  extent  where  analytical  predictions  of  the  residual  strength 
are  In  good  agreement  with  experimental  measurements  Q).  These  predictions 
are  based  on  transient  heat  transfer  analysis  of  the  thermal  heating  and 
straightforward  stress  analysis  to  determine  the  failure  load  from  a  fracture 
strength  versus  temperature  relationship.  In  the  Investigation  performed 
herein,  similar  analytical  methods  for  heat  flow  and  stress  analysis  have  been 
developed  to  calculate  the  crack  driving  force  In  flawed  metallic  structures. 
The  emphasis  of  the  methodology  will  be  to  establish  the  conditions  for 


nonduc+Ile  crack  extension  where  thermal  heating  Is  severe  enough  to  Increase 
tensile  loadings  In  critical  structure  components  where  cracking  may  exist. 

PROJECT  SCOPE 

In  failure  analysis  evaluations  of  flawed  structures,  the  application  of 
fracture  mechanics  provides  a  quantitative  means  of  assessing  structural 
Integrity.  The  fracture  mechanics  approach  to  structural  reliability  accepts 
that  flaws  will  exist  and  that  conditions  can  be  established  where  flaws  will 
remain  stable  and  not  grow  to  an  unacceptable  size  during  service  between 
Inspection  Intervals.  Fracture  mechanics  evaluations  require  the  calculation 
of  crack  tip  stress  Intensity  factor  (K)  which  defines  the  severity  of  the 
flaw  In  terms  of  Its  physical  size  and  the  applied  stress  acting  on  It.  By 
calculating  K  and  comparing  It  with  the  material  fracture  toughness,  the 
ability  of  a  given  loading  condition  and  flaw  dimensions  to  cause  unstable 
fracture  can  be  studied. 

The  development  of  a  general  methodology  for  determining  stress  Intensity 
factors  for  aircraft  structures  exposed  to  Intense  thermal  heating  will 
require  a  materials  property  data  base  at  elevated  temperatures,  a  nonlinear 
heat  transfer  analysis  model,  a  thermal  stress  analysis  model,  and  a  fracture 
mechanics  solution  technique  for  determining  K.  Although  the  problem  Is  very 
complex,  the  present  Investigation  analytically  examines  the  behavior  of 
through-cracked  and  part- through-cracked  aluminum  panels  subjected  to  Intense 
heating.  Several  simplifying  assumptions  and  model  Idealizations  have  been 
made  In  order  to  demonstrate  the  technique  In  this  Phase  I  project.  Stress 
Intensity  factors  are  computed  by  the  Influence  function  method  for  a  range  of 
crack  locations  within  the  panel  relative  to  the  position  of  the  heat  source. 
The  usefulness  of  the  method  Is  demonstrated  by  both  the  presentation  of 
numerical  results  and  the  Illustration  of  the  ability  to  solve  these  oompl  ex 
problems  on  desktop  microcomputers. 


Section  2 

STRATEGY  AND  OBJECTIVES 


The  strategy  employed  In  this  work  for  calculating  stress  Intensity  factors 
was  to  use  a  weight  or  Influence  function  approach.  As  will  be  discussed 
later,  this  technique  Is  very  efficient  and  enables  the  basic  method  to 
provide  accurate  results  using  small  microcomputers  without  the  need  for  large 
malnfrane  computers.  Because  the  Influence  function  method  requires  as  Input 
the  distribution  of  stress  In  the  structure  at  the  location  of  the  postulated 
flaw,  a  thermal  stress  solution  method  was  also  developed  which  can  be 
executed  on  a  small  microcomputer. 

The  primary  objectives  of  this  project  were: 

•  To  develop  an  Influence  function  algorithm  complete  with  the  thermal 
stress  solution  for  calculating  stress  Intensity  factors 

•  To  demonstrate  the  method  for  two  simple  engineering  problems 

•  To  specify  the  general  requirements  for  a  computer  program  that  will 
run  o"  a  desktop  microcomputer 

Feasibility  of  this  overall  approach  Is  shown  In  Section  5  by  successfully 
demonstrating  stress  Intensity  factor  calculations  for  the  following  two 
problems:  (1)  through-cracked  panel  subject  to  laser  beam  Impingement  In  the 
vicinity  of  the  crack  and  (2)  an  edge-cracked  pi  ate  subjected  to  aerodynamic 
heating  on  one  side.  Problem  Input  parameters  for  thermal  flux,  absorptivity, 
and  heat  losses  due  to  convection  and  radiation  were  thoroughly  researched  to 
provide  a  meaningful  demonstration  of  the  method.  Verification  of  the 
accuracy  of  the  developed  software  Is  shown  by  comparing  stress  Intensity 
factor  results  calculated  here  with  published  results  for  some  simple  cases. 


Section  3 

DETERMINATION  OF  STRESS  INTENSITY  FACTORS 


FRACTURE  MECHANICS  CONCEPTS 


itroductlon 


In  applying  fracture  mechanics  analysis  to  failure  prevention  evaluations.  It 
Is  Important  to  establish  the  possible  modes  by  which  the  structure  may  fall. 
Also,  the  parameters  which  are  Important  In  determining  the  residual  strength 
of  a  structure  containing  defects  must  be  defined.  The  failure  behavior  of 
structural  metals  can  be  classified  Into  three  regimes.  The  disciplines 
required  to  assess  these  regimes  are: 

•  Linear  elastic  fracture  mechanics  (LEFM)  -  The  structure  falls  In  a 
brittle  manner,  and  on  a  macroscale,  the  load  to  failure  occurs  within 
nominally  elastic  loading. 

•  Elastic-plastic  fracture  mechanics  (EPFM)  -  The  structure  falls  In  a 
ductile  manner,  and  significant  stable  crack  extension  by  tearing  may 
precede  ultimate  failure. 

•  Limit  load  or  plastic  collapse  -  The  failure  event  Is  characterized  by 
local  large  deflections  and  local  plastic  strains  associated  with 
ultimate  strength  collapse  at  a  cross  section  (the  structure  exhausts 
Its  redundancy  through  the  development  of  multiple  local  plastic 
Instabilities  until,  under  continued  application  of  load,  global 

col  I  apse  occurs) . 

A  schematic  diagram  showing  the  relationship  between  critical  or  failure 
stress  (i.e.  ,  residual  structural  strength)  and  flaw  size  Is  shown  In 
Figure  3-1  for  the  failure  modes  previously  described.  The  shape  and  position 
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Schematic  Showing  the  Relationship  Between  Failure  Stress 
and  Flaw  Size  For  Two  Limiting  Failure  Modes. 


of  the  failure  locus  will  depend  on  the  fracture  toughness  (K  )  and 

I  c 

strength  properties  (o  )  of  the  material  as  well  as  the  structural 

uts 

geometry  and  type  of  loading. 

Linear  elastic  fracture  mechanics  Is  used  most  appropriately  to  describe  the 
behavior  of  low  tough  ness/ high  strength  materials  In  which  the  plastic  zone  Is 
small  relative  to  the  structural  geometry  and  little  ductility  precedes 
fracture.  With  this  method,  no  account  Is  taken  of  Increased  material 
resistance  to  brittle  failure  when  significant  plasticity  occurs.  For  thermal 
stress  problems  considered  In  this  application,  It  Is  assumed  that  flaws  are 
primarily  driven  by  thermal  loads  that  will  elevate  stresses  local  to  the  flaw 
by  the  elastic  constraint  to  thermal  displacements.  Hence,  flaws  located  In 
the  vicinity  of  rapid  thermal  change  will  behave  elastically  provided  that  the 
thermal  stresses  are  nominally  low  (below  material  yield  strength).  For  this 
condition,  LEFM  theory  will  apply  and  we  restrict  our  discussion  to  this  case. 

Linear  Elastic  Fracture  Mechanics  Principles 

The  principles  of  LEFM  effectively  link  three  parameters — the  flaw  size,  the 
fracture  toughness  of  the  material,  and  the  applied  stress.  If  any  two  of 
these  are  known,  the  critical  value  of  the  third  can  be  quantified.  Although 
the  stress  distribution  of  a  cracked  structure  for  an  arbitrary  mode  of 
loading  and  shape  of  body  and  crack  can  be  quite  difficult  to  determine,  only 
three  deformation  modes  can  occur  near  the  tip  of  the  crack;  the  faces  can  be 
pulled  apart  (Mode  I)  or  sheared  perpendicular  or  parallel  to  the  leading  edge 
of  the  crack  (Modes  II  or  III).  These  three  loading  modes  are  shown 
schematically  In  Figure  3-2a  and  the  character  of  the  crack  near- tip  stress 
distribution  Is  Illustrated  In  Figure  3-2b.  The  crack  opening  mode  or  Mode  I 
In  which  the  load  Is  applied  normal  to  the  crack  face  Is  generally  considered 
the  most  damaging  of  the  three  modes. 

As  mentioned  above,  the  most  useful  parameter  for  describing  the  character  of 
the  crack  near-tlp  stress  distribution  Is  the  stress  Intensity  factor.  The 
stress  Intensity  factor  (K)  defines  the  local  crack  tip  response  to  global 
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(a)  Basic  Crack  Opening  Modes. 


All  stress  components  have  the  form: 

°ij  =  —  fii(k)(e) 

1J  vTivr  1J 

Where  i  r  x,  y,  z;  j  =  x,  y,  z;  and  k  =  1,  II,  III 


(b)  Near  Tip  Stress  Components. 


Figure  3-2  -  Review  of  Linear  Elastic  Fracture  Mechanics. 
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conditions  and  Is  calculated  In  terms  of  the  nominally  applied  stress  (o),  the 
crack  length  (a),  and  a  factor  that  depends  on  the  flaw  geometry,  stress 
distribution,  and  structural  displacement  constraints  (F(a))  from  the 
rel atlon: 

K  =  Fcv'na  (3-1) 

Assuming  Mode  I  loading,  fracture  Is  predicted  when  the  applied  value 

reaches  a  critical  level.  For  plane  strain  oondltlons,  this  critical  level  Is 

the  fracture  toughness,  K  ,  and  a  requirement  for  safe  service  Is: 

Ic 


K 
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Ic 


(3-2) 


The  critical  value  of  applied  stress  can  be  oomputed  In  terms  of  flaw  size  and 
fracture  toughness  from  the  expression: 


o 
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Fv-rra 


(3-3) 


Likewise,  given  the  applied  stress  and  critical  toughness,  the  critical  fl 
size  can  be  determined  Implicitly  from: 


a 


c 


(3-4) 


A  detailed  discussion  on  the  computation  of  K 


I 


fol  lows  next. 


ANALYSIS  fCTHODOLOGY 

Introduction 

Although  many  closed  form  or  approximate  solutions  exist  for  K,  often 
numerical  techniques  are  required  to  calculate  K  accurately  for  the  actual 
structure.  A  numerical  approach  was  required  for  the  thermal  problems 
considered  herein  because  of  the  complex  stress  state  expected  for  the  thermal 
behavior  of  metal  structures  subjected  to  rapid  thermal  heating.  Examination 


3-5 


of  Eq.  (3-1)  shows  that  the  key  element  In  calculating  K  Is  the  determination 
of  the  function  F(a),  which  carries  all  the  Information  concerning  the 
Influence  of  load  distribution  and  geometry.  Clearly,  If  K  has  been  computed 
as  a  function  of  o  and  a,  then  the  desired  nondlmenslonal  function  F(a)  Is 
determined  trivially. 

Tradltk  al  numerical  methods  for  calculating  K  are  either  energy  based  or 
crack  tip  stress/d  1 sp I acement  based  techniques.  These  stress  analysis  methods 
Involve  the  direct  modelling  of  crack  face  boundaries  with  very  refined 
regions  of  discrete  elements  or  nodes  In  the  stress  solution.  Such 
traditional  approaches  w 1 1 1  not  be  applicable  for  use  with  microcomputers 
because  of  the  significant  numerical  effort  required.  Under  certain 
circumstances,  solution  of  the  stress  Intensity  factor  for  one  set  of  loadings 
on  a  90I Id  provides  sufficient  Information  to  generate  easily  stress  Intensity 
factor  solutions  for  a  complete  class  of  loadings.  This  approach  Is  known  as 
the  Influence  function  method,  and  Its  highly  efficient  nature  Is  exploited  In 
this  Investigation.  This  characteristic  not  only  provides  accurate  K 
solutions  for  thermal  problems  but  also  allows  the  method  to  be  used  on  small 
computers. 

Influence  Function  Method 

The  Influence  function  method  Is  a  numerical  technique  that  allows  for  the 
calculation  of  K  for  nonlinear  (general)  varying  stress  distribution  acting  on 
the  crack.  The  Influence  function  or  weight  method  was  developed  by  Bueckner 
(3)  and  clarified  later  by  Rice  (4)  for  two-dimensional  problems.  The 
approach  was  expanded  to  three-dimensional  problems  by  Besuner  (5)  and  Cruse 
(fi).  Results  developed  by  Bueckner  In  1  958  ( J )  are  of  considerable  value  for 
application  to  problems  dealing  with  body  forces,  thermal  effects,  and 
residual  stress. 
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The  Influence  function  (h)  Is  a  function  of  crack  position  (x),  specified 
displacement  boundary  conditions  (u),  and  geometry.  The  calculation  of  K  for 
the  general  class  of  crack  problems  In  Mode  I  Is: 


-f 


h(x)o(x)dx 


(3-5) 


where  L  Is  the  crack  line  and  o(x)  Is  the  "uncracked"  stress  distribution 
normal  to  the  crack  face.  Once  the  Influence  function  has  been  formulated  for 
a  given  crack  configuration,  the  stress  Intensity  factor  for  any  applied 
stress  field  determined  from  the  uncracked  geometry  can  be  computed  by  simple 
numerical  Integration  of  Eq.  (3-5). 

The  essential  features  In  the  formulation  of  the  Influence  function  method  are 
based  on  the  following  fundamentals: 

e  The  application  of  elastic  superposition  allows  the  use  of  the 
uncracked  stress  distributions  In  the  K  analysis. 

e  The  Influence  function  Itself  Is  Invariant  with  stress  and  provides 
the  vehicle  to  calculate  the  effect  of  the  crack  In  red!  str  I  butl  ng  any 
stress  field. 

The  background  and  basis  of  the  above  principles  and  their  Impact  on  simply  Ing 
the  computational  effort  are  described  In  the  next  two  subsections. 


Superposition  and  Crack  Face  Loading  Equivalence 


The  principle  of  superposition  reduces  the  K  solution  of  an  arbitrary  and 
perhaps  difficult  crack  problem  to  two  simpler  problems— the  stress  analysis 
problem  but  without  the  crack  and  the  probl  an  of  a  cracked  body  with  an 
applied  pressure  that  cancels  the  uncracked  stress  field  to  establish  the 
traction  free  boundary  conditions  along  the  crack  face. 
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In  Bueckner's  formulation,  he  applies  the  theorems  of  Cl  apeyron  and  Betti  to 
demonstrate  that  K  for  an  arbitrary  cracked  body  subject  to  remotely  applied 
stresses,  displacements,  and  body  forces  Is  Identical  to  that  due  to  a  simple 
loading  system  on  the  crack  face  (2,  fi) •  This  principle  Is  Illustrated  In 
Figure  3-3.  The  general  cracked  body  problen  shown  In  Figure  3-3a  Is 
considered  to  be  the  sum  of  two  other  problems  shown  In  Figures  3-3b  and  3-3c. 


Hence: 


(b)  (c) 

K  +  K 


In  the  original  problem  (Figure  3— 3 a ) ,  the  cracked  body  Is  subject  to  loading 
Involving  surface  tractions  over  Boundary  S  ,  Imposed  displacements  over 
Boundary  S  ,  and  body  forces  within  the  volume  (V)  of  the  body.  A  stress 
free  crack  exists  In  the  body.  In  the  problem  given  In  Figure  3-3b,  the 
stress  distribution  prior  to  the  Introduction  of  the  crack  (the  "uncracked" 
stress)  Is  assumed  to  be  known,  and  the  surface  tractions  on  the  plane  where 
the  crack  Is  to  appear  are  denoted  as  T«.  The  tractions  T  *  are  stresses 
which  when  applied  to  the  crack  face  of  the  original  problem  are  Just 

sufficient  to  close  the  crack  completely.  Because  the  crack  Is  perfectly 

(b) 

closed,  K  equals  zero.  Hence,  the  model  problen  given  In  Figure  3-3c  Is 

Identical  to  that  of  the  original  problem  (Figure  3-3a).  In  this  case 
(Figure  3-3c),  the  body  Is  treated  with  zero  boundary  stresses  on  S  ,  fixed 
zero  displacements  on  S^,  no  body  forces  over  V,  and  only  acting  on 
the  crack  face.  The  strain  energy  (U)  Is  simply  the  work  done  by  T  *  on  the 
crack  face: 


■i  /  T'*u' 


(3-7) 


where  L  Is  the  crack  boundary.  The  rate  of  change  of  U  with  crack  length  for 
a  body  of  unit  thickness  Is  related  to  stress  Intensity  factor  by: 


1 U  _ 

3a 


(3-8) 
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The  parameter,  H,  represents  the  elastic  constants:  E  for  plane  stress  and 
2 

E/(1  -V  )  for  plane  strain  where  E  Is  modulus  of  elasticity  and  v  Is 
Poisson's  rat  to. 


As  will  be  shown  later  In  the  results  for  thermal  heating.  It  Is  highly 
desirable  to  view  the  original  crack  problem  of  Figure  3-3a  as  the  problem 
shown  In  Figure  3-3c.  It  Is  preferable  to  analyze  the  uncracked  body  to 
obtain  the  solution  for  T*  prior  to  solving  for  K.  This  uncracked  analysis 
permits  Identification  of  fl aw  locations  by  first  observing  high  stress  points 
In  the  body  as  well  as  determining  potential  crack  propagation  paths. 
Furthermore,  In  some  cases,  consideration  of  other  stresses  not  usually 
determined  by  numerical  stress  analysis  Is  required  (e.g.,  residual  stress). 
These  stresses,  however  determined,  can  be  simply  superimposed  and  substituted 
In  Eq.  (3-5). 


Influence  Function  Generation 


To  create  an  Influence  function,  h(x).  It  Is  necessary  to  obtain  the  solution 
for  T^«,  Uj,  and  K  for  a  simple  loading  condition  by  traditional  means  as 
suggested  by  Eqs.  (3-7)  and  (3-8)  and  to  differentiate  the  crack  face 
displacements  with  respect  to  crack  length.  Development  of  the  theoretical 
basis  Is  rather  complex  and  Involves  a  deep  understanding  of  analytic  function 
theories  (.5)  or  a  sophisticated  understanding  of  energy  principles  of 
elasticity  (4).  A  clear  depiction  of  the  theory  Is  possible  by  considering 
two  cases  of  loading  on  the  same  body  with  the  following  results: 


K  Is  the  stress  Intensity  factor  for  Case  1 


u  Is  the  displacement  vector  throughout  the  body  for  Case  1 


T  Is  the  stress  vector  appl  led  by  the  boundary  of  the  body  for 
Case  2 


(Z) 

•  F  Is  the  body  force  acting  within  the  body  for  Case  2 
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The  results  of  the  theoretical  analysis  (4)  show  that  for  Case  2  the  stress 
Intensity  factor  Is  given  by: 


2K  1  ^ 


.(2)  3u 


(2) 


(3-9) 


where  S  represents  the  boundary  of  the  body  and  A  Is  the  area  of  the  body. 
Since  the  solution  of  any  crack  problem  can  be  reduced  to  the  solution  of  the 
same  geometry  with  the  loading  on  the  crack  boundary  only,  we  can  confine 

_  _ ...  J  2)  _  (2)  (2).  . 


Eq.  (3-9)  to  the  case  where  F  =  0,  T 


equal  to  o  (x),  and  the 


only  displacement  of  relevance  In  Case  1  Is  the  displacement  normal  to  the 
(t) 

crack  face  u  over  the  region  of  the  crack  boundary.  The  expression 
(2)  n 

for  K  then  becomes: 


.(2)  . 
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(3-10) 


Comparison  of  Eq.  (3-10)  with  Eq.  (3-5)  gives  the  general  form  of  the 
Influence  function  as: 


h(x)  = 


“(TT  3  a” 


(3-11 ) 


Once  h(x),  the  Influence  function,  and  o(x),  the  "uncracked"  stress  acting 
along  the  hypothetical  crack  plane,  are  known,  the  stress  Intensity  factor  can 
be  calculated  by  simple  numerical  Integration.  This  Integration  can  easily  be 
accomplished  on  a  desktop  computer. 


INaUENCE  FUNCTIONS  FOR  TWO  CRACKED  BOOY  GEOMETRIES 


i 


Through  Cracked  Infinite  Plate 

The  Influence  function  for  a  through  crack  In  an  Infinite  plate  shown  In 
Figure  3-4  was  derived  from  the  uniform  tension  solutions  given  by  Irwin  (5) 
or  Erdogan  (IQ)  where  stress  functions  were  used  to  establish  the  following: 


K  -  a/na 


u(x) 


(2a  -  x)‘/2„ 


(3-1 2) 


From  Eq.  (3-11),  the  Influence  function  for  a  crack  over  the  region  0  <  x  <  2a 
under  Mode  I  loading  Is: 


for  a  crack  tip  located  at  x  =  2a,  and 

hj(x)  =  (na)"1/2  1/2  (3-1 3b) 

for  a  crack  tip  at  x  =  0.  The  Influence  function  for  Mode  II  and  Mode  III 
loadings  are  of  the  same  form  as  Eq.  (3-13),  hence: 


The  solution  for  K()  or  K  follows  from  Eq.  (3-5)  but  with  the  stress 

distributions  for  In-plane  (0  )  or  out-of-pl  ane  (0  )  shear  stress 

xz  xy 

distribution  substituted  for  the  normal  stress  0  (x). 


Figure  3-4  -  Through  Cracked  Infinite  Panel  of 
Uniform  Thickness. 


The  Influence  function  for  an  edge  cracked  plate  of  finite  thickness  was 
derived  by  Bueckner  (11)  for  the  geometry  shown  In  Figure  3-5.  A  series 
approximation  was  used  to  obtain  the  function  with  a  reported  accuracy  of 
1  percent  for  0  <  a/w  <0.5  (I.e. ,  crack  penetrations  up  to  one-half  the  plate 
thickness).  The  resulting  Influence  function  for  Mode  I  loading  Is: 


hr(x)  =  (2/t r)1/2  (a  -  x)1/2  1  +  Ml(l  -  -*)  +  M2(l  - 


(3-15) 


where, 

2  6 

M  =  A  +  B  (a/w)  +  C  (a/w) 

111  1 

2  6 

M  =  A  +  B  (a/w)  +  C  (a/w) 

2  2  2  2 

A  =  0.6147,  B  =  17.1844,  C  ■  8.7822 

1  1  1 

A^  =  0.2502,  Bz  =  3.2889,  C  «  70.0444 

For  the  case  of  edge  crack  In  a  half  space,  the  coefficients  become 

M  =  A  and  M  =  A  which  simplifies  the  expression  for  h  (x). 

112  2  I 

Equation  (3-15)  In  conjunction  with  the  general  stress  distribution,  o  (x), 

z 

can  be  Integrated  according  to  Eq.  (3-5)  for  0  <  x  <  a  to  obtain  the  stress 
Intensity  factor. 


Verification  of  Solutions 


The  solution  of  stress  Intensity  factor  for  the  two  Influence  functions  was 
accomplished  by  numerical  Integration  of  Eq.  (3-5).  A  simple  rectangular  rule 
was  used  to  perform  the  Integration.  Because  h(x)  Is  singular  at  x  =  a,  a 
nonuniform  Integration  grid  was  used  that  refines  the  Integration  steps  as  x 


Figure  3-5  -  Geometry  For  Edge  Cradled  Plate. 
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approaches  the  crack  tip.  Up  to  30  Integration  points  were  used  In  the 
numerical  procedure. 


The  accuracy  for  the  program  was  verified  by  comparing  the  numerical  results 
for  the  case  of  uniform  tension  to  that  of  known  values  given  In  the 
literature  (12).  This  comparison  Is  shown  In  Figure  3-6  for  both  flaw  models. 
Very  good  agreement  between  the  numerical  algorithm  and  the  literature  results 


I s  observed. 
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STRESS  INTENSITY  FACTOR,  F  = 


CRACK  LENGTH  OR  DEPTH,  a/b 


Figure  3-6  -  Comparison  of  Numerical  Results  From  Influence 
Function  Method  With  Theory  For  the  Case  of 
Uniform  Tension. 
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Section  4 
THERMAL  ANALYSIS  MODELS 


INTRODUCTION 

As  described  In  Section  3,  the  stress  Intensity  factor  model  requires  the 
thermal  stress  distribution  for  the  solid.  In  classical  thermal  stress 
analysis,  the  heat  transfer  solution  Is  decoupled  from  the  mechanics  solution 
since  thermal  deformations  usually  are  small  and  can  be  neglected  tn  the  heat 
transfer  analysis.  Although  In  severe  thermal  conditions  the  deformations 
could  be  large,  the  determination  of  stress  herein  will  be  restricted  to  small 
displacements  and  no  phase  change  so  that  the  decoupling  of  the  two  problems 
I  s  ecceptabl  e. 

The  classical  representation  of  heat  flow  In  a  three-dimensional  solid  Is 
given  by: 


(4-1) 


where  p  Is  the  material  density,  C  Is  the  specific  heat,  T  Is  the  temperature, 
t  Is  time,  and  k  Is  thermal  conductivity.  The  function  G  Is  the  energy 
produced  per  unit  volume  per  unit  time.  In  Eq.  (4-1),  p,  C,  and  k  are 
considered  functions  of  both  position  and  temperature  and  G  Is  a  function  of 
position  and  time.  This  equation  Is  a  statement  that  the  rate  at  which  energy 
accumulates  In  a  volume  Is  equal  to  the  net  flow  of  heat  across  the  surfaces 
of  that  vol  une  plus  the  rate  at  which  heat  Is  produced  within  the  volume. 

The  solution  of  Eq.  (4-1)  In  conjunction  with  the  appropriate  boundary 
conditions  and  Initial  conditions  gives  the  transient  temperature  response 
throughout  the  solid.  We,  therefore,  seek  solutions  of  Eq.  (4-1)  for  the 
cases  of  laser  Irradiation  and  aerodynamic  heating  so  that  the  desired  thermal 
stresses  can  be  oomputed. 
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One  of  the  most  Important  effects  of  Intense  heating  by  laser  Irradiation  Is 
the  conversion  of  the  electromagnetic  wave  energy  In  the  beam  Into  thermal 
energy  In  the  material  (JJ).  The  rate  and  source  of  production  of  heat  by  the 
laser  will  yield  the  function  G.  For  metals  exposed  to  a  continuous  beam  at  a 
wave  length  of  approximately  10  pm,  energy  absorption  occurs  In  a  very  thin 
layer  at  the  surface.  It  then  becomes  convenient  to  obtain  solutions  for 
Eq.  (4-1)  with  G  =  0  but  with  a  specified  heat  flux  at  the  surface  boundary. 
Hence,  for  continuous  wave  radiation  by  a  C0^  laser  (wave  length  of 
10.6  pm),  the  laser  beam  Impingement  problem  can  be  treated  similar  to 
aerodynamic  heating  In  that  the  Impinging  heat  flux  constitutes  a  boundary 
condition.  As  will  be  discussed  later,  other  surface  effects,  such  as 
reradlatlon  of  energy  and  heat  loss  due  to  convection  from  the  surface,  must 
also  be  Included  as  boundary  conditions. 


RADIAL  HEAT  TRANSFER  IN  A  SEMI- INFIN ITE  PLANE 
Governing  Equations 

The  representation  of  the  laser  beam  Impingement  plane  Is  modelled  as  a 
two-dimensional  radial  plane.  Figure  4-1  Illustrates  a  portion  of  a 
seml-Inf Inlte  metal  plate  of  thickness,  w,  with  Its  exposed  surface  Irradiated 
by  a  laser  beam  of  radius,  R.  A  cylindrical  coordinate  system  (r,  0  ,  z)  Is 
also  shown  with  the  z-axls  coincident  with  the  axis  of  the  laser  beam  and  the 
origin  located  at  the  metal  plate  back  wall. 

In  addition  to  the  Impinging  laser  beam,  the  exposed  surface  (z  =  w)  also 
experiences  two  other  modes  of  heating — convective  cooling  (or  heating)  by 
moving  air  at  all  values  of  r  and  surface  reradlatlon  at  all  values  of  r.  It 
Is  assumed  that  the  back  wall  (z  =  0)  Is  Insulated  at  all  values  of  r.  If  the 
laser  beam  Is  axl symnetrlc  about  z  and  the  oonvectlon  coefficient  Is  constant, 
then  the  resulting  heat  conduction  within  the  metal  Is  also  axl  symmetric  with 
the  result  that  the  metal  temperature  at  any  Instant  In  time  will  be  a 
function  of  r  and  z  but  not  angular  position  about  the  z  axis. 


W.VVVV  ’ 


CTC 


V  VV  V  * 

VO 


|>  V 


v;v;v>>.Vv>:\>v 


EXPOSED  SURFACE 


^  INSULATED  BACKWALL 

Metal  plate  thickness  w 

Radius  of  impinging  laser  beam  =  R 

Boundary  conditions: 

z  =  0,0<r<°°  +  Insulated 

z=w,  0  <  r  <  R  -►  Laser  beam  impingement 

0  s  r  <  00  -*■  Convective  cooling  and  surface  reradiation 


Figure  4-1  -  Schematic  of  Laser  Beam  Impingement  Area  on 
Exposed  Metal  Surface. 


The  time  varying  heat  transfer  In  the  metal  plate  Is  described  as  follows. 
Initially,  at  time  t  <  0,  the  plate  Is  at  a  uniform  temperature  throughout 
because  the  convective  cooling  (or  heating)  and  exposed  surface  rerad! atlon 
are  uniform  over  the  entire  surface.  Since  the  back  wall  Is  adiabatic,  there 
Is  no  localized  source  of  energy  anywhere  with  the  result  that  no  temperature 
gradients  are  promoted.  At  the  time  t  =  0,  the  laser  beam  strikes  the  exposed 
surface  over  the  circular  area  of  radius,  R.  The  beam  may  be  constant  or 
variable  with  time.  The  high  level  of  energy  associated  with  the  beam  Is 
conducted  down  Into  the  metal  and  outward  In  the  radial  direction.  Metal 
temperatures  at  z  =  w  and  r  <  R  rise  rapidly  Initially.  Temperatures  Indepth 
and  for  r  >  R  then  start  to  Increase  as  the  radial  conduction  commences. 
Counteracting  the  laser  beam  energy  source  are  the  two  energy  sinks  provided 
by  surface  reradlatlon  and  convective  ooollng. 

Due  to  the  heat  transfer  events  described  above,  the  temperature  field  In  the 
metal  plate  Is  two-dimensional  and  time  dependent,  T  =  T(r,  z,  t).  It  can  be 
shown  from  Eq.  (4-1)  that  this  temperature  field  In  cylindrical  coordinates 
with  no  Internal  heat  generation  Is  governed  by  the  following  partial 
differential  equation: 

_  3T  _  1  _3_  1T\  +  J  ft 

dC  3t  r  3r  Vk  3rJ  3t  3zJ  (4-2) 

where. 


p  =  p(T) 

=  Metal 

densl ty 

C  «  C(T) 

=  Metal 

spec! f Ic  heat 

k  =  k (T) 

=  Metal 

thermal  conductivity 

4-4 


Note  that  Eq.  (4-2)  Is  formulated  with  the  metal  thermophyslcal  properties  p, 
C,  and  k  allowed  to  vary  with  temperature.  Of  course,  Eq.  (4-2)  Is  subject  to 
both  an  Initial  condition: 


T  =  T(r,  z,  0) 


(4-3) 


for  all  z  and  r  values,  and  the  heat  flux  boundary  conditions  at  z  =  0  and 
z  =  w  already  described: 


T(r,  o,  t)  _ 
3z 


(4-4) 


.k  3T(r»  w>  t)  _  Q 
K  3z  g 


(4-5) 


where  0  Is  the  net  heat  flux  from  the  surface  energy  balance.  The  surface 
energy  balance  In  Eq.  (4-5)  Is  of  the  general  form: 

Q  =  q"(r,  w,  t)  =  -as  q  1a$er  +  q  c0nvection  +  q  radiation  (4-6) 

where  a  Is  the  surface  absorptance.  Thus,  for  the  laser  Impingement  region 
(0  <  r<  R): 


3T 


-k  v. 


3z 


-a  q\  +  h(T  -  T  )  +  a  ►  (T4  -  Tn4) 
s  H  laser  •*'  s  m  o 


(4-7) 


while  outside  of  the  laser  Impingement  region  (R  <  r  <  <»): 


3T 
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MT  -  u  *  os.  m(T4  -  T04: 


(4-8) 


where.  In  the  equations  above,  h  Is  the  surface  transfer  coefficient,  T«.  Is 

the  air  temperature,  T  is  the  remote  temperature  for  radiation,  o  Is  the 

o  s 

Stefan-Bol tzmann  Constant,  and  (  )s  5Urfoce  emlttance. 

ro 


Since  Eq.  (4-2)  Is  nonlinear  due  to  the  variation  of  thermophyslcal  properties 
with  temperature,  the  most  practical  means  of  solution  Is  to  use  a  finite 
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difference  method  which  Is  based  on  a  network  of  discrete  nodes  In  the  region 
of  the  metal  plate  where  significant  heat  conduction  occurs. 


Numerical  Solution 

In  this  work,  the  approach  described  by  Griffis,  et  al.,  (1_)  has  been 
followed.  Figure  4-2  Illustrates  a  lumped  mass  nodal  network  around  a  surface 
node.  This  node  experiences  heat  conduction  to  or  from  all  adjacent  Nodes  1, 
2,  and  3  and  heat  fluxes  to  or  from  Its  exposed  surface  due  to  the  three 
mechanisms  Identified.  As  derived  by  Griffis,  et  al . ,  (1),  the  Incremental 
temperature  rise  at  the  surface,  AT^,  for  an  Increment  In  time,  At,  Is  given 
by : 


1  *  Z  KiJ  (Tj  -  V  +  ri4riV"l.ser 


(4-9) 


*  ViVlA  '  Ti  1  +  ri4rih(T»  -  V  ri4r.A2'pjCi 


where,  for  the  radial  direction: 
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and,  for  the  vertical  direction: 


~Az.  Az  . 
k  '  k;’ 


(4-12) 


For  a  subsurface  node,  Eq.  (4-8)  simplifies  to: 


ATi 


LK.  .(T.  -  T.) 
1J  J  i 


L  J 


At 


r.Ar.AziPiCi 


(4-14) 


where  K  Is  defined  by  Eq.  (4-10)  for  the  radial  direction  and  by 
II 

Eq.  ( 4- 1  z )  for  the  vertical  direction.  Equations  (4-9)  and  (4-14)  are  the 
desired  finite  difference  approximations  to  Eq.  (4-2). 

As  discussed  by  Griffis,  et  al .  (1),  the  right  hand  side  of  Eq.  (4-9)  Is 
evaluated  at  present  time  and  the  left  hand  side  Is  treated  as  an  Incremental 
change  between  present  temperature  and  future  temperature  after  elapsed  time, 
At.  Mathematically,  this  Is  an  explicit  (rather  than  Implicit)  numerical 
solution  approach.  Because  of  the  thinness  of  aircraft  skin  relative  to  other 
dimensions,  the  through  thickness  heat  flow  Is  neglected  In  the  present 
Investigation.  This  assumption  will  allow  for  a  slmpl  Icatlon  of  the  method  to 
a  one-dlmenslonal  model. 


THROUGH  THICKNESS  HEAT  aCW  IN  A  CURVED  SHaL 

Governing  Equations 

The  thermal  response  resulting  from  one-dlmenslonal  heat  flow  through  the 
thickness  of  a  curved  cylindrical  shell  was  used  to  model  uniform  aerodynamic 
heating  of  the  leading  surface.  The  model  geometry  Is  shown  In  Figure  4-3. 
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(a)  Thermal  Model  Geometry  and  Boundary  Conditions. 
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(b)  Finite  Difference  Model. 


Figure  4-3  -  Thermal  Model  Geometry. 


The  solution  for  the  transient  response  of  two  concentric  shells  of  dlssmllar 
materials  Involves  solving  a  parabolic  partial  differential  equation  of  the 
form: 

92T  1  3T  _  J_  3T 

ar2  r  3r  at  3t  (4-15) 

where, 

t  =  Time 

r  =  Radial  coordinate 

T  =  Temperature,  T(r,  t) 

a  =  Thermal  dlffuslvlty 

t  7 

The  modelling  of  two-layer  construction  will  allow  the  analysis  of  aircraft 
structures  where  a  thermal  protective  skin  Is  appl  led  to  the  outer  surface. 

The  boundary  conditions  are  nonhomogeneous  and  require  that  shell  Inside 

(r  =  r  )  and  outside  (r  =  r  )  surfaces  transfer  heat  between  the  metal 
l  o 

surface  and  the  environment: 
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h  ,  h  =  Surface  heat  transfer  coefficients  (Inside  and  out- 
I  o 

side) 


k  ,  k  =  Thermal  conductivity  of  surface  materials 


T  ,  T  =  Metal  surface  temperatures 
wl  wo 


T  ^  =  Air  or  gas  temperature 


The  Initial  condition  for  the  analysis  Is: 


T(r,  0)  3  T  =  Constant  (4-18) 

0 


Numerical  Solution 

The  solutions  of  Eqs.  (4-15)  through  (4-18)  are  accomplished  using  a  finite 
difference  method  based  on  the  Crank-Nlcol son  technique.  This  Is  an  Implicit 
approach,  and  the  temperatures  are  obtained  by  solving  a  set  of  simultaneous 
linear  equations  (14).  The  nodal  point  scheme  for  the  difference  equations  Is 

A 

shown  In  Figure  4-3.  By  defining  the  temperature  T(r,  t)  =  T(R,  t)  -  T(r,  o), 

a  th 

the  difference  equations  T  will  represent  the  temperature  at  the  I 

I,  n  th 

point  Into  the  model  and  at  the  n  step  where  radial  distance,  r  =  lAr  and 
time  t  =  nAt.  For  simplicity,  the  WAn  notation  will  be  dropped  In  the 
difference  equations  which  will  follow. 

Two  additional  or  fictitious  node  points  are  used  so  that  a  central  difference 
relation  can  be  formulated  at  the  Inside  and  outside  shell  walls.  This  yields 
three  difference  equations  (one  for  each  node  at  the  wall  and  one  for  an 


Internal  node).  The  finite  difference  equation  for  the  convective  boundary 
condition  Is: 


+  1  ,n 

2Ar 


Ti,n> 


(4-1  9) 


The  finite  difference  equallon  for  the  partial  differential  equation  at  the 
wall  Including  the  convective  boundary  conditions  Is: 
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where  y(t)  Is: 
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The  finite  difference  equation  for  an  Internal  node,  say  the  I  +  3  node  as  an 
exampl  e.  Is: 
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For  the  case  when  an  Insulated  boundary  (h  =  05  Is  assumed,  the  finite 
difference  equation  for  the  Insulated  boundary  condition  Is  (central 
difference) s 


+  k  +  1  ,n  +  k  -  1  ,n  _  A 
2Ar  U 

This  results  In  a  difference  equation  for  the  node  at  the  wall  which  Includes 
the  Insulated  boundary  condition  as: 


Ti  +  k  -  1  ,n  +  1 


Ar 

atAt 


Ti  +  k,n+  1  +  Ti  +  k  -  1  ,n 


+  k ,n 


=  0 


(4-23) 


A  special  situation  arises  for  the  case  of  an  Internal  node  at  the  Interface 
of  fhe  dissimilar  materials.  For  this  node,  the  thermal  properties  of  the 
base  material  are  assumed  as  an  approximation,  although  a  more  exact  solution 
would  Involve  conservation  of  heat  transfer  across  the  boundary  between  the 
two  materials.  It  should  also  be  noted  that  the  temperatures  at  the 

convective  wall  (T  )  may  oscillate  In  time  about  a  central  ("correct")  value 

w 

when  the  heat  transfer  coefficient  Is  a  very  large  value.  This  phenomenon 
does  not  affect  the  accuracy  of  the  Interior  points  and  may  be  avoided  by 
using  a  smaller  time  Increment  (14). 


THERMAL  STRESS  ANALYSIS 


Once  the  temperature  distribution  Is  known,  the  elastic  stress  In  a 
two-dimensional  solid  can  be  determined  by  satisfying  the  equations  of 
equilibrium  and  compatibility  In  conjunction  with  the  boundary  conditions, 
treating  the  true  thermal  stress  problem  (l.e. ,  without  body  forces),  the 
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equlllbrlun  equations  will  automatically  be  satisfied  If  a  stress  function, 
can  be  found  such  that: 


d2^ 

~~Z  ’  °y 
dy 


-7»  ,  O 

2  xy 


(4-24) 


where  o  ,  o  ,  and  o  are  the  component  stresses, 
x  y  xy 

Considering  the  plane  stress  case,  the  strain  (€)  can  be  expressed  In  terms  of 
the  stress  function  through  the  use  of  stress  strain  relations.  Substitution 
of  these  strains  Into  the  compatibility  equation: 


32y 

_  _xy 


yields  the  Inhomogeneous  bl harmonic  equation: 


4  2 

V  *  +  Ea  VI  =  0 

e 


(4-25) 


(4-26) 


where  E  Is  the  modulus  of  elasticity,  a  fs  the  coefficient  of  thermal 
2  e 

expansion,  and  V  notation  represents  the  Laplaclan  operator.  Solution  of 
Eq.  (4-26)  for  the  stress  function  with  the  appropriate  boundary  conditions 
will  give  resultant  stresses  from  the  equlllbrlun  equations. 

Although  most  thermal  stress  problems  are  treated  by  approximate  solution 
methods,  a  class  of  exact  solutions  will  be  exploited  herein  to  simplify  the 
model  (15).  When  the  temperature  does  not  vary  In  one  direction.  It  can  be 
assumed  that  the  stress  due  to  heating  does  not  vary  In  that  direction.  For 
these  simple  cases.  It  Is  possible  to  Integrate  the  governing  stress  equations 
above  to  yield  the  thermal  stress.  For  the  circular  plate  geometry  for  laser 


values  for  T(r,  t)  computed  from  the  thermal  solution  routines.  For  fhe  model 
representing  aerodynamic  heating,  as  the  Interface  between  the  dissimilar 
materials,  another  node  Is  added  so  that  two  nodes  at  the  Interface  each  have 
the  material  properties  of  Its  respective  side.  Since,  In  general,  the 
temperature  distributions  across  the  Interface  will  be  continuous,  the 
accuracy  In  the  stress  calculation  at  the  Interface  will  not  be  In  significant 
error  even  though  a  considerable  stress  discontinuity  could  exist  at  the  same 


Section  5 

LASER  BEAM  IRRADIATION  ON  A  FLAT  PLANE 


PROBLEM  PARAMETERS 
Geometry 

The  problem  of  laser  beam  Impinging  perpendicular  on  a  flat  planar  surface  In 

the  vicinity  of  a  through  thickness  crack  Is  represented  In  Figure  5-1.  The 

crack  has  a  total  length  of  2a  and  the  center  of  the  crack  Is  a  distance, 

r  ,  from  the  center  of  the  laser  beam.  The  crack  Is  also  oriented  by  the 
c 

angle,  Q  ,  relative  to  a  radial  line  connecting  the  crack  center  with  the 
c 

laser  beam  Impingement  zone.  The  beam  cross  section  Is  assumed  circular  with 
a  radius,  R.  The  thickness  of  the  plane  Is  defined  as  w.  The  local 
coordinate  system  for  the  crack  Is  the  same  as  depicted  earlier  In  Figure  3-4. 

Because  of  the  many  variables  that  can  be  Investigated,  most  of  the  geometry 

parameters  are  held  fixed.  Specifically,  R  Is  equal  to  2  Inches,  w  Is  equal 

to  0.25  Inch,  and  the  ratio  r  /R  Is  set  at  two  In  the  analyses  that  follow. 

c 

The  effect  of  crack  length  and  Its  relative  angular  position  on  stress 
Intensity  factor  are  analyzed  parametrically. 


Thermal  Parameters  and  Material  Properties 

The  thermal  model  described  In  Section  4  requires  as  Input  the  surface  flux 
Imparted  by  the  laser,  the  absorptivity,  the  emlsslvlty,  and  the  basic  thermal 
properties  of  the  panel.  In  the  analysis,  the  panel  material  properties  are 
taken  from  typical  values  for  aluminum  alloy  Type  7075- T6  Q,  12,  IS).  A 
summary  of  the  assumed  properties  Is  given  In  Table  5-1.  The  thermal  and 
mechanical  properties  listed  are  at  room  temperature.  In  the  analysis,  the 
thermal  and  mechanical  properties  are  assumed  not  to  change  with  temperature. 
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Table  5-1 

SUMMARY  OF  PROPERTIES  OF  ALUMINUM  ALLOY  7075-T6 


Property 


Assumed  Value 


Specific  Heat,  C 
Density,  p 

Thermal  Conductivity,  k 
Thermal  Expansion,  ae 
Poisson's  Ratio,  v 
Elastic  Modulus,  E 
Yield  Strength,  oy 
Ultimate  Strength,  ou 
Liquidus  Temperature 
Solidus  Temperature 
Incipient  Melting  Temperature 


0.23  BTU/1  b-°F 
169  lbs/ft3 
75  BTU/hr-f t-°F 
13  x  10‘6  in/in-°F 
0.33 

10.3  x  106  psi 
73  ksi 
83  ksi 
1175°F 
890°F 
990°F 


fe*  i.<  ^  h*  t*i  fi  is 
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although  the  temperature  solution  algorithm  can  accommodate  temperature 
varying  properties. 

2  2 

Two  laser  radiation  fluxes  are  assumed  (100  W/cm  and  500  N/cm  )  which  are 

capable  of  being  achieved  by  a  large  C0^  laser.  The  surface  absorption  of 

the  aluminum  panel  can  range  between  0.6  to  0.9  depending  on  whether  the 

surface  Is  highly  polished  or  heavily  oxidized  or  painted  (13.  2Q).  In  this 

Investigation,  a  Is  assumed  to  be  0.8.  which  represents  a  painted  surface. 

Further  analysis  simplification  Is  made  by  assuming  gray  body  behavior  for 

radiation  so  that  emlsslvtty  and  absorptivity  can  be  equated.  Normally,  these 

parameters  are  surface  condition  and  wave  length  dependent  and.  In  general, 

Independent  of  each  other.  Remote  temperature  for  radiation  losses  from  the 

panel  Is  assumed  to  be  absolute  zero  (T  =  -460°F). 

o 


STRESS  SOLUTIONS 

The  temperature  and  stress  response  of  the  panel  was  solved  with  the  finite 
difference  model.  The  panel  was  modelled  In  one  dimension  with  only  a  single 
node  representing  the  global  z  (through  thickness)  direction.  A  total  of  200 
nodes  was  used  In  the  radial  direction  to  a  maximum  radius  of  20  Inches. 

Hence,  the  radial  spacing  between  nodes  was  0.1  Inch.  Because  of  the  Implicit 
numerical  procedure,  several  time  Increments  were  tried  In  order  to  verify 
numerical  stability.  It  was  determined  that  At  =  0.05  second  gave  reasonable 
results.  Transient  temperature  solutions  for  ten  seconds  duration  were  made, 
however,  because  of  the  analysis  assumptions,  only  the  solutions  up  to  four 
seconds  are  real Istjc. 

The  temperature  distribution  at  0.5  second  after  the  laser  flux  was  applied  Is 

shown  In  Figure  5-2.  In  this  example  and  the  calculations  for  K  that  follow. 

the  panel  was  assumed  stationary  so  heat  loss  due  to  convection  was  zero.  For 
2 

a  flux  of  500  N/cm  ,  the  metal  temperature  has  risen  by  more  than  200°F 
above  the  Initial  plate  temperature  of  70°F.  Because  the  laser  flux  Is 
applied  uniformly  over  the  region,  0  <  r/R  <  1,  the  temperature  profile 
exhibits  a  plateau  behavior  over  the  Irradiated  region  with  a  very  sharp 


•r.  -*^-r  w* 


attenuation  occurring  at  r/R  =  1.  The  case  of  100  M/cm  flux  shows  a 
significantly  less  severe  thermal  upset  of  the  plate  and  comparably  lower 
stresses  and  stress  Intensity  factor  would  result. 

The  radial  (o  )  and  tangential  (o^)  stresses  corresponding  to  the 
r  u 

temperature  gradients  of  Figure  5-2  are  plotted  In  Figure  5-3.  The  radial 

stresses  are  everywhere  compressive  Implying  that  cracks  oriented  parallel  to 

the  laser  bean  boundary  ((?  =  90°)  will  not  experience  positive  crack 

c 

opening  forces.  Based  on  the  magnitude  of  the  stress  components,  flaws 

oriented  perpendicular  to  the  beam  boundary  (.$  =  0)  will  have  the  highest 

c 

Mode  I  stress  Intensity  factors.  This  result  Is,  In  part,  due  to  our 
one-dlmenslonal  treatment  of  the  problem,  but  nonetheless,  this  observation 
will  be  true  In  most  two-dimensional  cases. 

The  highest  tensile  stresses  occur  close  but  Just  outside  the  beam  boundary 

and  attenuate  away  as  the  radial  distance  away  from  the  beam  Increases.  A 

time  equal  to  0.5  second,  the  peak  stress  Is  computed  to  be  9  ksl  for  the 
2  2 
500  W/cm  case  and  about  2  ksl  for  the  100  M/cm  flux.  Because  the  beam 

energy  Is  both  uniform  and  continuous  In  time,  the  stresses  continuously 

Increase  In  magnitude  under  the  elastic  conditions  assumed  In  the  model. 

STRESS  INTENSITY  FACTORS 

The  stress  Intensity  factor  was  calculated  by  the  Influence  function  method 
described  In  Section  3.  The  stress  Intensity  factor  for  a  crack  oriented 
perpendicular  to  the  beam  circular  boundary  Is  shown  In  Figure  5-4.  Here,  K 
Is  plotted  as  a  function  of  time  from  the  start  of  laser  beam  contact.  The 
crack  center  Is  located  4  Inches  from  center  to  the  beam  and  the  crack  half 
length  Is  1  Inch  (l.e. ,  a  =  1  Inch). 

The  K  level  at  each  crack  tip  Is  observed  to  Increase  with  time  as  would  be 
expected  for  a  continuously  local  heating  of  the  plate.  The  largest  K  values 
are  computed  for  the  highest  flux  case,  which  Is  also  expected.  What  Is 
Interesting,  however,  Is  that  the  K  level  due  to  heating  can  be  significant, 
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up  to  20  percent  of  K  for  the  material  and  occurring  within  the  first 

Ic 

second  of  Irradlatton  although  a  crack  of  2  Inches  In  total  length  would  be 
quite  large  to  remain  undetected  In  a  critical  airframe  member.  In  any  case, 
It  Is  possible  to  generate  large  crack  driving  forces  If  the  Impinging  beam  Is 
of  significant  energy  and  surface  area  contact. 

As  shown  In  Figure  5-4,  the  K  levels  for  each  crack  tip  are  very  nearly  equal 

early  for  the  thermal  transient  but  diverge  rapidly  as  time  progresses, 

2 

especially  for  the  500  M/cm  flux  case.  The  maximum  K  value  occurs  at  the 

crack  tip  (Crack  Tip  2)  closest  to  the  laser  heating  zone  due  to  the 

Increasing  stress  gradient  In  that  direction.  Unstable  fracture  Initiating  at 

that  crack  tip  would  undoubtly  propagate  toward  the  hot  spot  where  It  probably 

would  arrest.  Although  not  evaluated  In  this  Investigation,  the  conditions 

for  fracture  reinitiation  at  the  remote  crack  tip  can  also  be  determined  by 

selecting  appropriate  r  /R  ratios  (between  zero  and  two  In  this  example). 

c 

The  Influence  function  method  Is  capable  of  modelling  the  flaw  for  any  value 

of  r  . 
c 

The  stress  Intensity  factor  as  a  function  of  crack  length  for  the  stress 
conditions  occurring  at  2  seconds  for  the  two  laser  fluxes  Is  given  In 
Figure  5-5.  As  the  physical  crack  length  Is  Increased,  the  K  level  at  each 
crack  tip  also  Increases  with  the  crack  tip  closest  to  the  beam  exhibiting  the 
maximum  computed  K  value.  Although  not  calculated  In  this  Investigation,  It 
would  be  expected  that  flaws  located  closest  to  the  laser  Irradiated  zone  but 
not  directly  under  the  heated  region  would  be  the  most  critical  based  solely 
on  the  thermal  stress  distribution.  A  proper  structural  Integrity  evaluation 
would  combine  aerodynamic  loads  with  the  thermal  loads  In  order  to  establish 
the  critical  flaw  conditions.  This  feature  could  be  easily  added  to  the 
present  method. 


rr;.  For  the  case  of  a  crack  oriented  45°  with  respect  to  r  the  crack  tip 

deformation  mode  Is  predominately  Mode  II.  The  Mode  II  stress  Intensity 
«v 

<  factor  Is  plotted  as  a  function  of  crack  length  In  Figure  5-6.  Comparison  of 

f  % 

I  these  results  with  the  pure  Mode  I  case  of  Figure  5-5  shows  similar  trends  and 


these  results  with  the  pure  Mode  I  case  of  Figure  5-5  shows  similar  trends  and 

magnitudes.  For  the  45°  f  I  ww ,  K  Is  not  zero  but  Is  very  small  —  less  than 
1/2  1 

1  ksl  In  for  the  flaw  lengths  and  thermal  loads  shown  In  Figure  5-6. 

EFFECT  OF  AERODYNAMIC  COOLING 

In  the  previous  calculations,  the  flawed  panel  was  assumed  stationary  so  that 
the  surface  heat  transfer  coefficient  (h)  was  Input  as  zero  (no  forced  or 
natural  convection).  If  the  panel  was  moving  or  If  air  was  being  moved  over  a 
stationary  plate  at  the  given  velocity,  then  the  stresses  In  the  panel  would 
decrease  because  surface  ooollng  would  reduce  metal  expansion  In  the 
Irradiated  region.  Likewise,  the  stress  Intensity  factor  for  a  given  flaw 
condition  would  also  be  reduced. 

By  repeating  the  thermal  stress  analysis  with  convection  ooollng,  the 

2 

resulting  stress  Intensity  factor  for  a  flux  of  100  W/cm  Is  given  In 

Figure  5-7.  In  the  thermal  analysis,  the  heat  transfer  coefficient  was 

estimated  frcm  the  case  of  uniform  air  flow  over  a  flat  plate  assuming  a 

turbulent  boundary  layer.  At  a  free  stream  velocity  of  1,100  ft/s  (Mach  1), 

2 

the  heat  transfer  coefficient  was  computed  to  be  115  BTU/hr-ft  -°F.  In 
Figure  5-7,  the  ratio  of  K  to  the  K  level  computed  for  zero  velocity  (h  =  0) 

Is  shown  as  a  function  of  velocities  up  to  Mach  1.  For  Irradiation  exposure 
of  1  second,  the  reduction  In  K  Is  approximately  4  percent  at  Mach  1.  These 
reduction  ratios  are  approximately  the  same  for  all  flaw  lengths  due  to  the 
short  exposure  time.  After  4  seconds,  however,  K  Is  reduced  more  with  a 
greater  reduction  observed  for  longer  flaws  (a  =  1  Inch)  than  shorter  flaws 
(a  =  0.25  Inch)  as  shown  In  Figure  5-6.  The  maximum  difference  between  K  for 
stationary  plate  to  K  I  n  a  plate  moving  at  Mach  1  Is  about  10  percent  for  the 
conditions  depicted  In  Figure  5-6.  It  Is  expected  that  the  benefits  of 
aerodynamic  cooling  will  become  more  significant  for  longer  exposures  and 
higher  fluxes  where  temperatures  closer  to  melting  could  occur.  Clearly,  the 
maximum  excursion  In  K  would  be  limited  by  the  presence  of  surface  ooollng  for 
the  case  of  a  single  laser  pulse  since  the  total  heat  absorbed  would  be 
reduced  by  convection  as  opposed  to  Just  heat  conduction  alone. 


STRESS  INTENSITY  FACTOR,  K..  (KSI  - 1 N  */2  ) 


Figure  5-6  -  Mode  II  Stress  Intensity  Factor  for  Two  Laser 
Flux  Cases  After  2  Seconds  of  Exposure 
(Note:  K.  <  1  k s i  in^  for  the  Above  45 
Crack ) . 
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Section  6 

AERODYNAMIC  HEATING  OF  A  PLANE  SURFACE 


PROBLEM  PARAMETERS 

A  second  problem  Involving  the  surface  heating  of  a  finite  thick  plate  was 
solved  In  similar  fashion  to  the  laser  beam  problem.  In  this  example,  the 
edge  cracked  plate  model  was  used  to  Investigate  the  effect  through  thickness 
temperature  gradients  have  on  a  two-dimensional  surface  flaw.  The  thermal 
model  Illustrated  In  Figure  4-3  was  used  to  9olve  for  temperatures  and 
stresses  as  a  function  of  time.  The  plate  was  assumed  to  be  one  material  (no 
cladding)  with  material  properties  the  same  as  given  In  Table  5-1  and  with  a 
thickness  of  1  Inch. 

Aerodynamic  heating  Is  applied  uniformly  to  one  side  (r  *  0)  and  the  flaw  Is 
located  on  the  other  side  (r  *  w)  where  the  surface  Is  assumed  adiabatic.  The 
surface  heat  flux  was  modelled  as  a  traditional  convective  heating  with 
surface  reradiation  effects  associated  with  a  black  body  radiator.  An 
equivalent  heat  transfer  coefficient  and  convective  gas  temperature  were 
estimated  from  measurements  on  metallic  specimens  subjected  to  space  shuttle 


reentry  conditions  (21).  Maximum  heat  fluxes  of  about  50  BTU/ft  -sec 
2 

(57  W/cm  )  and  a  surface  heat  transfer  coefficient 
-3  2 

0.008  x  10  Ibm/ft  -sec  based  on  the  enthalpy  gradient  between  the  wall 
and  the  air.  An  equivalent  h  and  T  were  established  from  these  data  by 
equating  surface  heat  flux  with  estimated  air  temperatures  at  corresponding 
free  stream  test  pressures  and  enthalpies  to  give: 


h  =  18  BUI/ ( hr- ft  °F 


T®  =  10,000°F 


(6-1 ) 


y*  ,-v  ,v  jv.-.-.vw.v. 
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It  should  be  noted  that  the  resulting  heat  rate  Is  approximately  an  order  of 
magnitude  less  severe  than  the  laser  heating  problem  analyzed  In  Section  5. 


STRESS  SOLUTIONS 

The  temperature  and  stress  solution  for  the  plate  was  solved  by  the  Implicit 
finite  difference  model.  The  plate  was  represented  by  21  nodes  with  a  linear 
spacing  of  0.05  Inch.  A  constant  time  Increment  of  0.1  second  was  used  In  all 
transient  solutions.  The  temperature  change  for  times  up  to  30  seconds  Is 
shown  In  Figure  6-1.  The  Initial  temperature  of  the  plate  was  70°F.  At 
1  second,  the  surface  temperature  was  computed  to  be  approximately  130°F,  but 
the  temperature  at  x  =  w  still  remains  unelevated.  At  30  seconds,  the  AT 
across  the  plate  Is  approximately  100°F. 

The  tangential  (o  )  stress  for  the  transient  times  of  0.1,  1,  and  10  seconds 
z 

Is  shown  In  Figure  6-2.  Although  the  stress  variation  with  distance  through 
the  thickness  Is  nonlinear,  the  basic  response  of  the  plate  will  be  to  bend 
like  a  beam  due  to  the  thermal  upset  at  the  surface.  The  through  thickness 
stress  Is  negative  and  relatively  small  compared  to  o  and,  therefore.  It  Is 
not  plotted.  Inspection  of  the  stress  distributions  to  30  seconds  Indicates 
that  the  maximum  stress  In  the  region  where  the  flaw  will  be  postulated 
(x/w  =  1)  occurs  between  7  and  10  seconds.  The  peak  tensile  surface  stress 
occurs  at  9.4  seconds. 


STRESS  INTENSITY  FACTORS 

Again,  the  Influence  function  method  was  used  to  calculate  the  Mode  I  stress 
Intensity  factor.  The  K  variation  as  a  function  of  time  Is  given  In 
Figure  6-3.  As  suggested  by  the  stress  solutions,  the  peak  stress  Intensity 
factor  occurs  at  approximately  9  seconds  after  heating  begins  and  slowing 
attenuates  with  time.  Ftve  K  curves  are  plotted,  each  corresponding  to  a 
different  but  constant  crack  penetration  (e/w  =  0.1  through  0.5).  Maximum  K 

level  for  each  flaw  depth  Is  observed  to  occur  at  approximately  the  same  point 

1/2 

In  time  with  the  largest  value  being  16.9  ksl-ln  for  a  flaw  that  Is 


HEATED 
SIDE  " 


DISTANCE,  Vw 


Figure  6-2  -  Thermal  Stress  in  a  Slab  Subjected  to  Aerodynamic 
Heating  on  One  Side. 
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Figure  6-3  - 


Stress  Intensity  Factor  Variation  with  Time. 
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Section  7 

COMPUTER  HARDWARE  AND  SOFTWARE  REQUIREMENTS 


INTRODUCTION 

The  oomputer  programs  developed  for  this  Investlgtton  were  all  written  In 
FORTRAN  and  were  executed  on  an  IBM-3081  system  located  at  Stanford 
University.  The  FORTRAN  language  was  selected  because  of  Its  general  use  In 
engineering.  These  programs  were  written  to  demonstrate  the  Influence 
function  method  and  to  obtain  approximate  estimates  of  execution  times  for 
various  machines.  These  oomputer  codes  are  not  considered  to  be  production 
programs.  Further  developments  will  be  required  before  the  basic  algorithms 
are  suitable  for  general  use.  Recommendations  for  additional  capabilities  are 
given  In  Section  8. 

The  purpose  of  this  section  Is  to  discuss  the  applicability  of  these  methods 
for  use  on  desktop  computers.  Since  the  major  concern  here  will  be  accuracy 
and  speed,  these  two  basic  characteristics  were  studied.  In  addition,  a 
specification  for  a  general  purpose  program  Is  given. 


PROBLEM  EXECUTION 

The  problems  of  laser  Impingement  and  aerodynamic  heating  were  solved  In  two 
basic  steps.  First,  the  thermal  stress  solution  (heat  flow  and  stress 
analyses)  was  performed  followed  by  the  stress  Intensity  factor  solution. 
Each  step  was  performed  by  an  Individual  program.  The  performance  of 
different  computers  was  determined  by  comparing  the  execution  times  of 
Individual  microcomputers  with  the  central  processing  unit  (CPU)  execution 
time  of  the  IBM-3081  system. 

The  microcomputer  systems  studied  were  all  based  on  Intel  8088  or  80286 


microprocessors  In  a  system  architecture  similar  to  IBM-PC/XT  and  IBM-PC/AT 
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systems.  Because  clock  speed  and  the  use  of  a  math  coprocessor  can  greatly 
enhance  machine  performance,  these  hardware  features  were  also  tested.  The 
coprocessor  for  the  80  88  CPU  Is  the  8087  and  the  coprocessor  used  with  the 
80286  CPU  Is  the  80287.  All  machines  contained  and  were  executed  from  hard 
dl sks. 

A  summary  of  Individual  machine  performance  was  as  measured  by  execution  time 
as  shown  In  Table  7-1.  All  times  listed  are  In  units  of  seconds.  Three 

problems  were  solved  as  shown — a  thermal  stress  solution  for  laser  heating,  a 

thermal  stress  solution  for  aerodynamic  heating,  and  a  stress  Intensity  factor 
for  an  edge  cracked  plate.  Execution  times  Included  both  CPU  time  and 
Input/output  times  associated  with  the  disk.  Although  the  slowest  algorithm 
appears  to  be  the  explicit  finite  difference  solution  for  thermal  stress,  this 
problem  required  many  more  nodes  In  the  model.  Hence,  comparison  between 

algorithms  should  not  be  made  In  Table  7-1.  It  would  be  expected  that  the  run 

times  for  the  Implicit  method  would  significantly  Increase  with  the  number  of 
nodes  used;  however,  the  stability  of  the  method  would  allow  for  larger  time 
Increments  to  be  employed. 

All  execution  times  displayed  by  the  microcomputers  are  within  the  realm  of 
practical  application;  however,  waiting  up  to  35  minutes  for  a  thermal 
solution  and  20  minutes  for  corresponding  K  solution  for  a  8088  machine 
without  a  coprocessor  will  undoubtedly  be  annoying  to  the  user.  The 
microcomputers  exhibited  execution  times  that  ranged  from  40  to  1,400  times 
slower  than  the  IBM-3081  mainframe  system.  The  greatest  Improvement  Is 
observed  when  a  coprocessor  Is  used  where  a  factor  of  five  on  speed 
enhancement  was  recorded.  This  Is  consistent  with  the  math  Intensive 
numerical  procedure  and  double  precision  arithmetic  used  In  the  algorithms.  A 
40  percent  speed  Improvement  was  observed  when  clock  speed  was  changed  from 
4.77  Wz  to  8  MHz. 

As  expected,  the  80286  machine  was  the  best  performing  system  with  run  times 
less  than  1  minute.  This  class  of  machines  Is  well  suited  for  the 
applications  Intended. 
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f-dl  stress  solutions  were  executed  with  the  WATF1V  Compiler.  The  K  solutions  were  executed  on  FORTRAN  Level  G  Compiler, 
rt  thermal  analysis  Tor  200  nodes  for  100  time  increments  with  stresses  calculated  for  16  time  increments, 
nt  thermal  ana’ysis  for  21  nodes  for  1 00  time  increments  with  stresses  calculated  for  16  time  increments, 
t.  i  or.  of  seven  '  versus  crark  depth  distributions  each  containing  50  discrete  values  for  t. . 


V. 


The  accuracy  of  microcomputers  as  compared  with  the  results  from  the  IBM-3081 
was  excellent.  The  word  length  for  the  8088  system  was  determined  to  be  24 
bits  for  the  mantissa.  Therefore,  the  8088  microcomputer  has  the  same  number 
of  significant  digits  In  floating  point  arithmetic  as  the  mainframe.  In  fact, 
the  mantissa  for  single  precision  arithmetic  In  the  3081  system  ranged  between 
21  bits  to  24  bits  depending  on  the  size  of  the  number.  This  suggests  that 
the  microcomputers  studied  herein,  on  average,  have  more  precision  than  the 
3081  mainframe. 


SOFTWARE  SPECIFICATION 

The  software  Is  not  organized  for  production  work.  Improvements  will  be 
required  In  the  area  of  Input/output  formats,  solution  efficiency,  and  post 
processing  of  results.  A  flowchart  showing  the  general  numerical  procedure  Is 
shown  In  Figure  7-1.  It  Is  proposed  that  FORTRAN  be  used  wherever  possible 
because  of  Its  general  use  as  a  scientific  programming  language.  The  Input  to 
the  progran  will  require  the  greatest  development  because  all  well  accepted 
programs  are  those  that  are  simple  to  use,  to  understand,  and  to  modify  Input 
files.  Programs  that  can  accept  both  batch  input  files  as  well  as  menu  driven 
Input  have  the  greatest  flexibility  to  create  new  problem  Input  files  or 
change  existing  or  previously  created  Input  files.  The  use  of  "help  files" 
will  greatly  Improve  the  speed  of  solving  problems  when  a  users  manual  Is  not 
conveniently  available  to  answer  questions. 


j 
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The  solution  algorithms  can  follow  the  basic  structure  of  the  programs. 

Because  these  programs  were  written  for  solving  demonstration  problems,  their 
efficiency  Is  not  optimum.  Additional  software  will  be  required  to  provide 
for  a  library  of  properties  for  appropriate  aircraft  structural  materials. 

The  properties  should  cover  the  temperature  range  expected  for  the  application 
of  the  method.  The  Influence  function  flaw  models  should  also  be  organized  In 
a  library  format  for  ease  of  selection  and  use  by  the  analyst. 

Recommendations  on  expanding  the  capabilities  of  the  present  demonstration 
software  are  discussed  In  Section  8. 
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The  processing  of  output  wilt  take  the  form  of  either  printed  summary  tables 
or  graphs  as  selected  by  the  user.  Either  output  format  option  should  have 
the  capability  of  being  displayed  on  the  monitor  as  well  as  printed  or  plotted 
on  a  hard  copy  device.  Displaying  of  results  on  the  oomputer  monitor  Is  very 
useful  In  quick  problem  solving  and  Input  Iterations,  especially  If  they  can 
be  performed  In  a  tlmeshare  environment.  Both  tlmeshare  and  batch 
environments  are  easy  to  Implement  on  a  single-user  microcomputer  system. 


Section  8 

SUWARY  AND  CONCLUSIONS 


The  Influence  function  method  Is  a  very  powerful  technique  for  calculating 
stress  Intensity  factors  given  that  the  stress  state  for  the  structure  without 
the  flaw  Is  known.  For  thermal  stress  problems,  the  Influence  function  method 
Is  very  useful  because  repetitive  K  solutions  are  often  required  due  to  time 
varying  stress  conditions  and  multiple  flaw  orientations  of  Interest.  While  a 
knowledge  of  the  stresses  along  postulated  crack  planes  Is  required  for  Input 
to  the  method,  these  stress  solutions  can  be  solved  numerically  by  algorithms 
that  complement  Influence  function  methods  In  both  simplicity  and  speed.  This 
capability  has  been  demonstrated  for  two  Intense  heating  problems — laser  beam 
Impingement  and  aerodynamic  heating.  Results  from  these  analyses  suggest  that 
stress  Intensity  factors  for  the  thermal  conditions  assumed  can  be  significant 
and  must  be  accounted  for  In  structural  reliability  assessments  where  severe 


surface  heating  Is  an  expected  or  potential  event. 


The  computer  resources  required  to  perform  these  computations  were  not 
significant,  and  all  computations  can  be  performed  on  currently  available 
desktop  microcomputers.  Although  execution  times  for  microcomputers  were  very 
slow  compared  to  an  IBM-3081  mainframe  system  (up  to  1,400  times  slower 
depending  on  microcomputer  architecture),  problem  run  times  of  less  than 
2  minutes  are  expected  for  a  80286-based  microcomputer,  such  as  an  IBM-PC/AT 
personal  computer.  Therefore,  the  Influence  function  method.  Including 
supporting  thermal  heat  flow  and  stress  analysis,  can  be  accommodated  by 
currently  available  microcomputers  that  already  enjoy  wide  acceptance  and  use. 
Also,  the  accuracy  of  the  methods  did  not  degrade  when  executed  on  the  smaller 
microcomputers.  The  methods  developed  herein  are,  therefore,  candidates  for 
use  as  a  desktop  design  tool  as  well  as  an  Inf  I Ight  expert  system  for  real 
time  damage  assessment. 
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Because  the  software  developed  during  this  project  was  for  the  purpose  of 
demonstrating  the  methods,  the  programs  are  not  set  up  for  production  problem 
execution.  Additional  refinements  to  the  methods  are  required  before  the 
methods  can  be  used  on  general  thermal  problems.  Specifically,  the  present 
programs  should  be  enhanced  In  the  following  areas: 

•  Expand  the  thermal  heat  flow  and  thermal  stress  numerical  procedure  to 
two-dimensional  analysis 

•  Improve  the  thermal  heat  flow  algorithm  for  speed  and  stability; 
explore  more  stable  methods  that  are  explicit  In  nature 

•  Include  temperature  varying  properties  and  surface  phase  change 
(melting,  sublimation) 

•  Expand  the  nunber  of  Influence  functions  to  Include  flaw  models  for 
cracks  emanating  from  structural  discontinuities,  such  as  holes  and 
stiffener  attachments 

•  Provide  for  ooupling  between  thermal  solution  and  stress  solution  when 
large  deflections  occur 

•  Verify  accuracy  of  overall  thermal  analysis  procedure  by  performing 
analytical  and  experimental  test  problems 

•  Explore  the  application  of  the  computer  software  and  hardware  deve¬ 
loped  here  to  an  aircraft  onboard  expert  system 

•  Prepare  user  documentation  for  the  method 

The  laser  Impingement  problem  Illustrated  that  local  plasticity  within  the 
radiation  zone  and  vaporization  are  possible  In  extreme  cases.  At  a  mlnlmun, 
these  conditions  should  be  checked  by  the  program  and  the  solution  procedure 
modified  to  account  for  this.  These  advanced  problems  can  be  handled  by  the 


method  provided  that  the  plasticity  Is  contained  In  the  zone  and  the  flaw  Is 
well  surrounded  by  an  elastic  medlun.  For  the  case  of  complete  burn  through, 
the  situation  becomes  the  problem  of  a  crack  in  the  vicinity  of  an  expanding 
hole  that  Is  being  heated  on  Its  boundary.  Such  a  case  oould  still  be  solved 
by  the  methods  described  herein,  provided  that  the  proper  boundary  value 
problem  solution  for  stress  Is  developed. 

In  conclusion,  the  objectives  of  the  project  have  been  achieved.  The 
Influence  function  method  Is  a  viable  technique  for  determining  stress 
Intensity  factors  for  realistic  thermal  boundary  assumptions.  Because  of  fhe 
simple  numerical  procedure,  the  method  will  execute  efficiently  on  small 
microcomputers. 
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